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ABSTRACT 

cn 

We propose a model for how the buildup of dark halos by merging satellites 
(Nj 1 produces a characteristic inner cusp, of a density profile p oc r~ ain with a in — > 

a as £ 1, as seen in cosmological N-body simulations of hierarchical clustering 
scenarios. Dekel, Devor & Hetzroni (2003) argue that a flat core of < 1 exerts 
tidal compression which prevents local deposit of satellite material; the satellite 



(N 



sinks intact into the halo center thus causing a rapid steepening to a- in > 1. Using 



merger N-body simulations, we learn that this cusp is stable under a sequence 
of mergers, and derive a practical tidal mass-transfer recipe in regions where the 
local slope of the halo profile is a > 1. According to this recipe, the ratio of mean 
densities of halo and initial satellite within the tidal radius equals a given function 
which is significantly smaller than unity (compared to being ~ 1 according 
to crude resonance criteria) and is a decreasing function of a. This decrease 
& ■ makes the tidal mass transfer relatively more efficient at larger a, which means 

steepening when a is small and flattening when a is large, thus causing converges 
to a stable solution. Given this mass-transfer recipe, linear perturbation analysis, 

> ■ 

supported by toy simulations, shows that a sequence of cosmological mergers with 
homologous satellites slowly leads to a fixed-point cusp with an asymptotic slope 
a as > 1. The slope depends only weakly on the fluctuation power spectrum, 
in agreement with cosmological simulations. During a long interim period the 
profile has an NFW-like shape, with a cusp of 1 < < a as . Thus, a cusp is 
enforced if enough compact satellite remnants make it intact into the inner halo. 
In order to maintain a flat core, satellites must be disrupted outside the core, 
possibly as a result of a modest puffing up due to baryonic feedback. 

Subject headings: cosmology: theory — dark matter — galaxies: formation — 
galaxies: halos — galaxies: interactions — galaxies: kinematics and dynamics 



1. INTRODUCTION 



A relatively robust universal shape for the density profile of dark-matter halos has 
been seen in cosmological N-body simulations of dissipationless hierarchical clustering from 
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Gaussian initial fluctuations. It can be approximated by the functional form 



P(r) = Ps 



( 
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where r s and p s are characteristic inner radius and density respectively. 1 This density profile 
is characterized by an inner "cusp" oc r~ a[n and a continuous steepening through a bend 
near r s towards r -3 near the virial radius R v (defined by a fixed mean overdensity A v above 
the universal mean, with A v = 180 to 340, depending on time and the cosmological model). 
Navarro, Frenk & White (1995; 1996; 1997, hereafter NFW) found eq. (1) with a in ~ 1 and 
a out ~ 3 to be a good fit to halos in simulations over the radius range (0.01 — l)R v , for a wide 
range of halo masses and for a range of hierarchical cosmological scenarios with different 
power spectra of initial fluctuations. Cole & Lacey (1996) came to a similar conclusion 
for self-similar scenarios with power-law power spectra, Pk oc k n with n = 0,-1,-2, in 
an Einstein-deSitter cosmology. High-resolution simulations of a few individual halos in a 
cosmological environment (Moore et al. 1998; Ghigna et al. 2000; Klypin et al. 2001) found 
that the typical asymptotic cusp profile at r < r s is sometimes somewhat steeper, closer 
to a in ~ 1.5. A careful convergence analysis by Power et al. (2002), who explored the 
robustness to numerical errors, found for the standard ACDM cosmology that a in reaches 
a slope shallower than 1.2 at their innermost resolved point of r ~ 0.005i? v . Thus, the 
robust result of the simulations is that dark halos have inner cusps with a characteristic 
slope 1 < a- m < 1.5. We seek a basic theoretical understanding of the origin of these cusps. 

An even more intriguing puzzle is introduced by observations which indicate that at 
least in some cases the actual inner halo density profiles are close to flat cores, with a m ~ 0. 
This has been detected directly by rotation curves in low surface brightness (LSB) galaxies, 
whose centers are dominated by their dark halos (van den Bosch et al. 2000; de Block et al. 
2001; Marchesini et al. 2002). Cores have also been argued to exist in normal disk galaxies, 
where more involved modeling is required (Salucci & Burkert 2000; Salucci 2001; Borriello 
& Salucci 2001). The presence of a core seems to introduce a severe challenge to the CDM 
cosmological paradigm. We mention in the conclusion section several mechanisms that have 
been proposed for the origin of such cores. In particular, attempts to turn a cusp into a 
core by direct stellar feedback effects in the present halos, which looked promising at a first 
sight (Navarro, Eke & Frenk 1996), seem not to work (e.g., Geyer & Burkert 2001; Gnedin 
& Zhao 2002). 

We find it useful to first try to understand in simple basic terms the possible origin of the 
universal cusp in the gravitational N-body simulations of cold dark matter. This may provide 



a useful generalized functional form has been proposed by Zhao (1996) 
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us with a tool for addressing the formation and survival of flat cores by other mechanisms, 
in particular by baryonic feedback processes within the hierarchical CDM framework. 

Several mechanisms have been studied in the context of dark-halo profiles in dissipation- 
less simulations. The outer slope of r~ 3 (and steeper) may possibly be explained in terms of 
violent relaxation (e.g., Barnes & Hernquist 1991; Pearce, Thomas & Couchman 1993 and 
references therein). We note in general that any finite system would tend to have a steep 
density fall off at large radii due to diffusion of particles outwards. Secondary spherical infall 
is expected to produce a profile closer to p oc r~ 2 , which may explain the behavior in the 
intermediate regions of the halo, but is too steep to explain the flatter inner cusp (Lokas & 
Hoffman 2000 and references therein). Thus, none of the above mechanisms seem to provide 
a natural explanation for the characteristic cusp of a ^ 1. 

The clear impression from the cosmological N-body simulations of hierarchical clustering 
scenarios is that halos are largely built up by a sequence of mergers of smaller structures. 2 
In a typical merger, a bound satellite halo spirals into the center of the larger halo due to 
gravity and dynamical friction. The satellite gradually transfers mass into the host halo due 
to tidal stripping or by eventually melting into the halo inner region. This process is likely 
to have an important effect in shaping up the density profile. Indeed, Syer & White (1998), 
Nusser & Sheth (1999) and Subramanian, Cen & Ostriker (2000, hereafter SCO) argued, 
using certain simple models and simulations of a sequence of mergers, that the buildup by 
mergers may naturally lead to a stable profile. However, they find their predicted profile to 
be quite sensitive to the power spectrum of fluctuations and to allow an inner slope of a < 1, 
in conflict with the robust result of the cosmological simulations. In fact, when trying to 
repeat the Syer & White analysis using their simplified modeling of the stripping process but 
with higher resolution, we find that in the long run the profile does not really converge to a 
stable cusp but rather continues to steepen slowly towards a m = 3 (see §8). Either way, it 
seems that something is not adequate in the simplified model adopted to describe the mass 
transfer from the satellite to the halo. 

We re- visit the buildup of halo profile by merging satellites and gain an encouraging 
new insight. We add two important new ingredients to the tidal effects. In another paper 
(Dekel, Devor & Hetzroni 2003, hereafter DDH) we argue that for a flat mean density profile 
with a < 1 the tidal effects on typical satellites induce three-dimensional compression with 
no local mass deposit, which results in a rapid steepening of the inner profile to a > 1. In 
this paper, we derive a useful prescription for tidal mass transfer at a > 1, and obtain higher 
deposit efficiency at higher a. We then show that this tends to flatten steep profiles with 



2 though cusps may also emerge in other scenarions where mergers are suppressed, see a discussion in §8. 
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large a and thus slowly leads to an asymptotic fixed point at a certain a = a as £ 1. 

In §2 we address the compression at a < 1 and the resultant steepening (DDH). In §3 
we derive a simple mass-transfer prescription in the range a > 1, based on merger N-body 
simulations and toy-model understanding. In §4 we explain why this prescription should lead 
to an asymptotic cusp. In §5 we use linear perturbation analysis to compute the asymptotic 
slope for satellites of a given mass, and in §6 we extend the analysis to a cosmological 
distribution of masses. In §7 we demonstrate this process via a semi-analytic simulation of 
a cosmological sequence of mergers. In §8 we discuss our results. 



This introductory section describes the robust transition of a core (a in < 1) to a cusp 
(ctin > 1), which is the main theme of another paper, DDH. We summarize it here for 
completeness and as a background to the independent analysis described in the following 
sections. 

A useful quantity in describing the tidal forces exerted by a halo of mass profile M{r) 
is the local logarithmic slope of its mean density profile p(r) oc M(r)/r 3 , 



such that locally p oc r~ a . We assume that a is either constant or monotonically increasing 
as a function of r, with values in the range < a < 3. The extreme values of a = and 3 
correspond to a const ant- density halo and a point mass respectively. Note that if the profile 
inside r is a power law, then the local and mean density profiles have the same logarithmic 
slope. In general, they are related via p{r) = [1 — a(r)/3]p(r), but they do not necessarily 
have the same slope at a given r. The slope of p(r) is equal to or larger than the slope of 
p(r). The following analysis refers to a as the slope of p(r). 

We then consider a satellite of mass m v <C M v , moving under the gravity excreted by 
the halo, when its center of mass is at position r as measured from the halo center. The tidal 
acceleration exerted by the halo mass distribution on a satellite particle at position vector £ 
relative to the satellite center of mass is obtained by transforming the gravitational attraction 
exerted by the halo on the particle into the (non-rotating) rest frame of the accelerated 
satellite. In the tidal limit £ <C r, to first order in £/r, using Cartesian coordinates about 
the satellite center, where £ = (£i, £2, £3) and £\ lies along r, we obtain 



2. 



CORE TO CUSP BY TIDAL COMPRESSION 




(3) 
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The components perpendicular to the line connecting the centers of mass are always of 
compression towards the satellite center, with an amplitude that does not explicitly depend 
on a. The maximum radial tidal force outwards is obtained along the line connecting the 
centers of mass. In the limit where the tides are exerted by a point-mass halo, a — 3, the 
pull outwards is maximal. For flatter halo slopes, the tidal stretching becomes weaker in 
proportion to (a — 1), until it vanishes at a — 1 and reverses direction into compression for 
a < 1. Thus, while for a > 1 there is always a tidal component pulling outwards, for a < 1 
the tidal forces are of compression everywhere in the satellite, and in the limit of a core with 
a = the tides induce symmetric compression in all directions. 

We argue in DDH that this critical transition at a ~ 1 is the main source for the origin 
of a cusp steeper than r _1 . The idea is that if the local tidal mass transfer from the satellite 
to the halo stops when the satellite's orbit has decayed into a core region where ct(r) < 1, 
the satellite would continue to sink in due to dynamical friction without further mass loss 
until it settles in the halo center. This would inevitably cause a general steepening of the 
core profile towards a > 1. In DDH we show in different ways that, indeed, no tidal transfer 
of mass from the satellite to the halo is expected in a region where a < 1. We address this 
point in simple cases using analytic approximations in the impulse and adiabatic limits and 
then demonstrate the anticipated effects using merger N-body simulations. 



Once satellites continue to merge with a halo of an inner cusp a > 1, the final slope is 
determined by the details of the tidal mass transfer from the satellite to the halo. We analyze 
this process via a quite general toy model prescription, which we justify and calibrate using 
N-body simulations. 



We consider a halo of mean density profile p{r) oc M(r)/r 3 with a corresponding log- 
arithmic slope profile a(r), and a merging satellite of initial mass profile m{£) with a cor- 
responding mean density profile a(£) oc m(£)/£ 3 . Let m f (r) describe the final distribution 
of stripped satellite mass in spheres about the halo center. By equating m f (r) and m(£) we 
obtain a one to one correspondence between £ and r. We then define the initial-density ratio 



3. MASS-TRANSFER PRESCRIPTION AT a > 1 



3.1. Toy Model Prescription 




(4) 
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This is an operational definition independent of any model assumption. We argue below that 
ip(a) is a relatively robust function quite insensitive to the specific nature of the merger, with 
the distinct property that it is a monotonically decreasing function of a, and with values 
significantly below unity in the range a > 1. The robustness of ip(a) is indicated in §3.2 
by a preliminary set of merger simulations, to be confirmed and refined by a more complete 
suite of simulations in an associated paper. 

In order to obtain a qualitative feeling for why ip{ot) should decrease with a we appeal 
to a very simplified toy model. We assume that the halo is spherical and fixed during the 
merger, and that the satellite is spherical and is tidally stripped shell by shell outside a 
momentary tidal radius. The stripped material from a given satellite shell, identified by its 
original radius £, is assumed to be deposited on average in a halo radius r. 

The condition ip — 1 refers to a crude resonance condition, a(£) = p(r), where the 
orbital period at £ within the (unperturbed) satellite equals the orbital period of the satellite 
within the halo at r. This is commonly assumed to approximate the momentary tidal radius 
(e.g., Syer & White 1998; Klypin 1999a). However, this condition ignores the inevitable a 
dependence of the stripping process, and it fails to address the structural changes of the 
satellite before stripping and the difference between where the stripping occurs and where 
the mass is actually deposited. 

We demonstrated in DDH that the tidal mass transfer becomes weak at a < 1, which 
means that ip is higher there compared to where a > 1. If we define the tidal radius £ t by 
the Lagrangian point where the net force in the satellite rest frame vanishes along the line 
connecting the centers of mass of halo and satellite, we obtain for a satellite on a circular 
orbit 



If the satellite structure inside the tidal radius is assumed to remain fixed, and if the stripped 
material outside £ t is assumed to be effectively deposited at the r where it is stripped, we 
obtain ip(a) oc a^ 1 and oc (a — l)^ 1 respectively. This is the kind of decrease in ip(a) 
expected due to the a dependence of the stripping efficiency 

Another source of a dependence is due to the systematic difference between the radii 
of stripping and deposit. Along a typical orbit, the satellite distance from the halo cen- 
ter oscillates periodically between the radii of apocenter r a and pericenter r p , while their 
amplitudes gradually decay due to dynamical friction. Ghigna et al. (1998) studied the 
distribution of satellite orbits in a high-resolution N-body simulation of a cluster emerging 




(5) 



and for a satellite on a radial orbit 



[a(r) - 1] p(r) = a(£ t ) . 



(6) 
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from a CDM cosmological background and found that the median ratio r a /r p is 6:1, with 
radial orbits common and circular orbits rare, and with a distribution of eccentricities quite 
independent of r. They also demonstrated that the tidal radii of the satellites are consistent 
with being determined near pericenter, approximately under the general resonance condi- 
tion p(r p ) ~ er(-^t)- Particles that escape from the satellite can be assumed, on average, to 
continue on an orbit about the halo center with apocenter and pericenter radii "frozen" at 
their values near the time of escape, suffering no further decay due to dynamical friction. 
Since a particle spends most of its time near the apocenter of its orbit, we can say that the 
escapers are effectively "deposited" near the apocenter radius valid at the time of escape. 
We can thus assume that the ratio of stripping radius ( £ r p ) and deposit radius ( £ r a ) 
is typically e £ r p /r a ~ 1/6. If we ignore the structural evolution of the satellite inside 
the tidal radius as well as the a dependence of the stripping efficiency addressed above, we 
obtain straightforwardly a decrease as a function of a of the sort ip(a) = e a , where a is the 
effective slope of the mean density profile between r p and r a . 

A third source of a dependence is the distortion of the satellite inside the tidal radius, 
with an a-dependent stretching along one axis and an a-independent compression along 
the others. The resultant decrease in satellite density before stripping has been seen in 
simulations (e.g. Klypin et al. 1999a, Fig. 6; Hayashi et al. 2002). At higher a values the 
stretching is stronger, which should make the stripping even more efficient there. 

The above three effects provide qualitative hints for the expected properties of ip(a), 
and in particular for it being a decreasing function of a with values below unity. We turn 
to N-body simulations in order to quantify these predictions. 

3.2. Merger Simulations 

In order to evaluate ip(cx) at a > 1 and test its robust properties, we ran several N- 
body simulations of isolated mergers between a large halo and a satellite halo of mass ratio 
m/M = 0.1. Similar simulations were used in DDH to study the inner core region and 
they are also described there using several additional figures which are not essential for our 
current purpose. We use the Tree code by Mihos & Hernquist (1996 and references therein) 
but with dark- matter halos only (no gaseous disks). In our default simulations, the large 
host halo is represented by 10 5 equal-mass particles and the satellite by 10 4 particles. We 
re-ran one case with 0.55 x 10 6 particles and the force resolution higher by a factor v^5, 
and confirmed that the results were practically identical. The simulation units are: length 
3.5kpc, mass 5.6 x 1O 1O M , and time 13.06Myr. The force softening length is 0.08 units, i.e. 
0.28kpc 
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The initial halo density profile, as measured in the unperturbed initial conditions, is a 
truncated isothermal sphere with a flat core, 

with p s = 10.36, r s = 1 and r\ = 10 The density at the characteristic radius r s is thus 
p(r s ) = 5.13. The internal velocities are constructed to fulfill the isotropic Jeans equation 
which ensures an equilibrium configuration as discussed in Mihos & Hernquist (1996). When 
run in isolation, the halo profile has been tested to be very stable for many dynamical 
times. The logarithmic slope a(r) of the mean density profile spans the range of interest 
between a = and 3, and its variation as a function of radius can be described to a good 
approximation by a(r) « 1.73 logr + 0.67 throughout the range 0.3 < a < 2.9. The initial 
halo density profile can be seen in Fig. 3 below. It resembles the generalized NFW profile of 
eq. (1) with a core of a <C 1. We chose to start with an inner core in order to also simulate 
the rapid steepening to a cusp, but the outer profile, where a > 1, is quite generic. 

The satellite initial density profile is fit by a Hernquist profile, 

with the default choice a s = 19.2 and 4 = 1 defining our typical "compact" satellite. The 
initial satellite density profile can also be seen in Fig. 3. In the inner region it convergence to 
the NFW profile, a- in =1. If we fit this Hernquist profile with an NFW profile by matching 
the characteristic radii where the local logarithmic slope is —2, we find that the radius 
corresponding to the characteristic radius of the NFW profile is I = 4/2. At this radius the 
density is 16o~ s /27, which for our default satellite is a(£ s /2) = 11.38. 

The satellite parameters were chosen to roughly mimic a typical compact satellite ac- 
cording to the distribution of halo properties in the ACDM scenario. In a hierarchical 
clustering scenario with a fluctuation power spectrum Pk oc k n , the halo characteristic radii 
and densities scale like £ s /r s oc m( 1+ ^/ 3 and cr s / p s oc mT v ', with v — (3 + n)/2. For ACDM 
on galactic scales one has v ~ 0.33. This is confirmed by cosmological simulations (NFW; 
Bullock et al. 2001a), where the typical halo profile is NFW. The halos of lower masses are 
thus typically more compact; they tend to have lower characteristic radii and higher cor- 
responding mean densities within these radii (corresponding to higher virial concentration 
parameters). For CDM halos of mass ratio m/M = 1/10, the typical ratio of characteristic 
radii and corresponding densities are expected to be roughly 0.4 and 2.1. The corresponding 
ratios in the initial conditions of our simulations are approximately 0.5 and 2.2, providing a 
reasonable match. It may also be interesting to note that the satellite mean density interior 
to 4/2 (or 4) is about 3.9 (or 1.09) times the halo mean density interior to r s . 
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Fig. 1. — The 10 4 satellite particles in 6 snapshots during the elongated merger, projected 
onto the orbital plane. The 10 5 live halo particles are not shown. The center of mass is at 
the origin and the dot marks the temporary halo maximum density. The circle is of radius 
r = 20 ~ 2r t about the halo maximum density, corresponding to where the initial halo 
practically ends and where the satellite is at the onset of the simulation. The snapshots 
correspond to the initial conditions, the first pericenter, the following three apocenters and 
the final distribution of stripped satellite material. The three thirds of the mass, in concentric 
shells about the satellite bound center at each time, are marked by different colors. The 
distance of the satellite center from the halo center is marked by r (while radii within the 
satellite are marked by £). 

We simulated three cases of initial merger orbits: a radial orbit, a circular orbit, and a 
typical elongated orbit with r p /r a ~ 1/6 initially. The unperturbed satellite is put initially 
at r = 20, i.e. at about 2r t , where the initial circular period is about 230Myr. In the circular 
and radial cases the magnitude of the initial satellite velocity was set to equal the circular 
velocity of the halo at that radius, namely a bound orbit with an orbital kinetic energy 
that equals half the absolute value of the total energy. For the elongated orbit the initial 
tangential velocity was one half of the circular velocity at r = 20. Each merger has been 
followed until the satellite's bound core has practically settled at the halo center. 

Fig. 1 shows the satellite mass distribution in 6 snapshots during the elongated merger, 
projected onto the orbital plane (while DDH show analogous plots for the radial and circular 
mergers). The satellite oscillates about the halo center through a sequence of pericenter and 
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apocenter passages. The shape of the orbit remains roughly constant while it shrinks in 
scale; when measuring the ratio of pericenter to the following apocenter it is r p /r a ~ 1/3.5, 
and when measuring the ratio of pericenter to the preceding apocenter it is r p /r a ~ 1/3.5, 
but these ratios remain roughly the same for all detectable pericenters. The oscillations 
decay due to dynamical friction until the satellite becomes confined to the halo core after 
~ 125Myr and 5 pericenter passages. By the first pericenter, the satellite is already stretched 
and stripped along its orbit, while it is temporarily shrunk in the perpendicular direction (as 
expected in DDH, see §2). This is followed by a re-bounce and significant mass loss about 
the following apocenter. The visual impression confirms the notion that the particles that 
are torn away near a pericenter radius continue on orbits that reflect on average the satellite 
orbit at the time of stripping, while the bound remnant continues to sink into smaller radii 
due to dynamical friction. For example, in the snapshots corresponding to the second and 
third apocenters we clearly see a large amount of stripped satellite material spread about the 
location of the previous apocenter. We can crudely say that mass is stripped near pericenter 
and is practically "deposited" about the following apocenter radius. The final distribution 
of satellite mass extends quite smoothly about the halo center in a puffy oblate ellipsoid that 
looks quite symmetric in the orbital plane. 

Fig. 2 describes the time evolution of the satellite spherical mass profile during the 
elongated merger by showing the mean radii of concentric spherical shells about the satellite 
maximum-density center, each encompassing a given fraction of the satellite mass (i.e., not 
necessarily the same population of particles at different times). Pericenter passages are 
identified five times. Overall contraction of bound satellite shells seems to start roughly at 
these times, as expected both in the impulse and adiabatic limits (see DDH). Each major 
contraction is followed by a re-bounce as the satellite moves towards apocenter, which results 
in overall expansion and stripping of the outer shells. This seems like a manifestation of 
the expected delayed stripping due to the energy pumped into the satellite by the inwards 
impulse in the core, as well as the adiabatic stretching and stripping outside the tidal radius. 
Once the satellite becomes confined to the inner halo where a < 1, after ~ 125Myr and 5 
apocenters, there is no apparent overall shell expansion anymore, indicating that stripping 
has stopped, as expected. Similar effects are noticed in DDH for the radial and circular 
mergers. 

Fig. 3 addresses the cusp formation straightforwardly by showing the density profile of 
the halo before and after the merger on a radial orbit. In DDH we also show similar plots 
for the circular and elongated merger orbits, with very similar results. For the standard, 
compact satellite, the left panel demonstrates the inevitable steepening of the profile in the 
core region, starting near the core boundary where a ~ 1. The slight depletion of final 
density in the region near a ~ 1, compared to the slight increment in halo density at larger 
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Fig. 2. — Time evolution of the satellite spherical mass profile during the elongated merger. 
Shown in the bottom panels are the mean radii of concentric spherical shells about the 
momentary satellite center, each encompassing a given fraction of the original satellite mass, 
as marked on the right. The top panels show the time evolution of position r of the satellite 
center relative to the halo center, and the corresponding local slope a(r). The times of 
pericenter passage are marked in the bottom panel by vertical lines. The stripping of a shell 
can be crudely identified by a rapid increase of its radius. Lack of stripping (and even a 
slight overall contraction of the innermost bound shells) is noticed whenever the satellite 
enters the halo core. 



radii, is consistent with no mass transfer in this region while the orbit of the remaining 
satellite continues to decay into smaller radii. 

How stable is this cusp under further mergers with similar satellites? To test this, we 
performed a series of mergers following each other. The main progenitor halo in each merger 
is taken to be the outcome of the previous merger, except that its total mass is scaled down 
back to the original halo mass (such that all the density profiles can be directly compared). 
We achieve this by letting each simulation involve the same numbers of N halo particles 
and n satellite particles (each of a fixed mass m), where the halo particles are selected at 
random from the N + n particles of the halo produced in the previous merger. In these 
tandem mergers we used N = 20, 000 and n — 2, 000. The halo profile in the range of 
interest is found to remain stable for at least several hundred Myr when the halo is run in 
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Fig. 3. — Halo density profile before (solid) and after (dashed) the simulated mergers with 
satellites of mass ratio 1:10. The initial halo has a flat core. Left: one merger with a 
compact satellite. A significant fraction of the satellite settles intact at the halo center 
without depositing mass near a £ 1, causing the steepening of the core into a cusp. Middle: 
tandem mergers with compact satellites on random elongated orbits. After each merger, the 
total mass of the halo is scaled down back to the original mass. The cusp remains stable 
under the subsequent mergers. Right: one merger with a puffy satellite, whose density has 
been scaled down by a factor of 3.3 compared to the compact satellite. The puffy satellite 
loses most of its mass in the outer halo, leaving the halo core practically unaffected. 

isolation. The satellite starts with the same profile as the original compact satellite. It is 
put at the same distance of r = 20 as before but in an elongated orbit of a random spatial 
orientation. The initial satellite velocity is tangential to the line connecting the centers of 
mass of halo and satellite, with an amplitude in units of the circular velocity at r = 20 chosen 
at random in the range v/v c = 0.2 — 0.6 (compared to the typical case simulated before, 
where v/v c = 0.49 led to peri/apocenter ratio of roughly 1:6). Each merger was followed 
for 325Myr before the following merger started. The density profiles after each of the first 
5 mergers are shown in the middle panel of Fig. 3. The first merger reproduces a cusp very 
similar to the cusp produced when the merger was simulated with N = 100, 000 (and with 
N = 500, 000), indicating that even N = 20, 000 is adequate for crude results in the regime 
of interest here. We see that inside the initial core radius the cusp is stable at a slope of 
a ~ 1.5, with no systematic tendency to deviate from the power-law profile resulting already 
after the first merger. At the radius where a ^ 1.5 the profile flattens, and the amplitude 
becomes lower partly because of the renormalization of the total mass between every two 
mergers. Our main goal in the current paper is to try to understand the origin of the stable 
cusp as indicated by these simulations of tandem mergers. 

How much puffing is needed in order to prevent damage by merging satellites to the 
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halo core? In order to obtain a first clue, we have performed simulations similar to those 
described above, except that the compact satellites of CDM have been replaced by more 
puffy satellites of the same total mass. In the Hernquist profile, eq. (8), where the default 
compact satellites had a s = 19.2 and £ s = 1, the parameters of the puffy satellites are now 
<7g = <7 s /8 and £' s = 2£ s . With this choice, the density at £ = 1/2, the characteristic radius 
of the original compact satellite, is scaled down by a factor of 2.3, to just below the density 
of the host halo at its characteristic radius. This corresponds to a reduction by a factor of 
2.8 in the mean density interior to £ = 1/2. The right panel of Fig. 3 shows the effect of 
such a merger with a puffy satellite on the halo density profile in the case of a radial radial 
merger (DDH also show similar results for other merger orbits). When the satellite is puffy, 
we see that almost all the satellite mass is stripped before the satellite orbit decays to the 
a < 1 zone, and as a result the halo core is practically unaffected. We learn that a modest 
reduction in the initial satellite inner density is enough for preventing the cusp formation 
seen in the compact-satellite case. 



3.3. Measuring ip(a) in the Simulations 

The measurement of ip(a) in the simulations is straightforward (and free of any model 
assumption). We simply measure the mass profile of the stripped satellite mass about the 
halo center at the final time, mf(r), and equate it with the initial satellite mass profile m(£) 
to obtain the deposit relation £(r). Then ip is evaluated from eq. (4) at any desired r, and 
is expressed as ip(a) given the halo slope profile a(r). 

Fig. 4 shows the measured ip(a) in three different merger simulations, where the orbits 
are very different from each other: circular, radial and elongated 1:6. All the curves show 
a similar general behavior, consistent with our qualitative expectations. First, the value of 
ip is significantly below unity, demonstrating that the mass transfer is more efficient than 
predicted by the naive ip = 1 model. Second, the function ip(a) is monotonically decreasing 
with increasing a in the range 1 < a < 3. The general behavior can be crudely approximated 
by several different functions, such as the power law ip(a) ~ 0.5/a, or the exponential 
i/j(a) ~ 0.5°. The scatter about these fits is roughly ±0.1. The relative robustness of ip(a) 
to the nature of the merger orbit indicates that eq. (4) may serve as a useful approximate 
recipe for tidal mass transfer in a general merger. 

In order to evaluate the sensitivity of our results to the resolution of the simulation, 
we ran a case identical to the merger on an elongated orbit described above, but now with 
5 times more particles (N = 0.55 x 10 6 ) and a softening length smaller by a factor of 5 1 / 3 
accordingly. The results are found to be practically identical in all respects, namely the decay 
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Fig. 4. — The mass-transfer condition. Density ratio ip, at deposit, as a function of halo 
slope a. Shown are the measured results from the three different N-body merger simulations. 
A qualitative general fit in the range a > 1 is provided by the power law ip(a) ~ 0.5/a, or 
by the exponential i/j(a) ~ 0.5°. 

rate of the satellite radius within the halo, the mass loss from the satellite, the final halo 
density profile, and the derived ip{a). This indicates that the resolution of our simulations 
is adequate in the range of radii of relevance, well inside the core/cusp and down to below 
0.2r s . 

In the current simulations we have only tested in a limited way the relative robustness to 
the merger orbit. The simulated halo and satellite are of a typical mass ratio and their profiles 
relate to each other in general accordance with the expected average scaling in the ACDM 
cosmology, which makes the obtained ^(a) a sensible first guess. However, the robustness of 
i/j(a) to variations about this mass ratio and these profiles should be confirmed and refined 
using a more complete suite of merger simulations. We adopt below ip(a) = 0.5/a as an 
illustrative example for obtaining numerical results but the following analysis is valid for a 
quite general ip(a). 

4. AN ASYMPTOTIC PROFILE: QUALITATIVE 

We saw in DDH that when the halo profile is flat, a < 1, tidal compression causes rapid 
steepening to a > 1. In this section we investigate the development of the profile due to a 
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sequence of mergers between similar halos where a is of order unity or larger. We show that 
if the tidal stripping is described by a condition similar to eq. (4), with ip(a) monotonically 
decreasing rapidly enough, then the profile evolves slowly towards an asymptotic stable 
power law r _Q?as , with a as larger than unity, which we can evaluate once ip(a) is given. The 
evolution to this asymptotic profile is through a sequence of profiles which crudely resemble 
the generalized NFW shape, and with an inner cusp of slope slightly larger than unity, 
1 < a < a as . While the value of a as may depend on the exact shape of ip(a), the following 
analysis showing the convergence to an asymptotic profile is not sensitive to it as long as 
i]){a) is decreasing rapidly enough. 



4.1. Homologous Halos 

We consider the buildup of the halo of mass M by a sequence of mergers with satellites 
of masses m < M drawn from a cosmological distribution of halos (similar to Syre & White 
1998). In our search for a self-similar evolution, we assume that the halo and satellite before 
the encounter are homologous, meaning that their unperturbed mean density profiles are 
scaled version of each other, 

<t(e) = 6 p(e/\), Q = m/\\ (9) 

where m stands for the satellite-halo mass ratio m/M (we use hereafter M — 1). If we 
refer, for example, to the functional form of eq. (1), then g = o~ s / p s and A = £ s /r s . The 
cosmological N-body simulations of CDM as well as of power-law power spectra (NFW; 
Bullock et al. 2001a) show that halos of lower masses tend to have lower characteristic radii 
and higher corresponding mean densities within these radii (corresponding to higher virial 
concentration parameters). The scaling of the averages of these quantities, at any given 
time, is found to be 

A = m {1+u)/3 , Q = rrT v , (10) 

with v ~ 0.33 for ACDM. Given a power-law power spectrum P(k) oc k n , a simple scaling 
argument based on linear theory predicts this behavior with v — (3 + n)/2. Thus, v ~ 0.33 
corresponds ton~ —2.3, as expected for a ACDM power spectrum at the relevant scales. 



4.2. Convergence to an Asymptotic Profile 

The origin of an asymptotic slope can be qualitatively understood in simple terms 
via the illustration in Fig. 5. We show a schematic mean density profile of a halo and a 
homologous satellite scaled accordingly. As in the toy- model interpretation of eq. (4), we 
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Fig. 5. — A schematic illustration of satellite mass deposit in the halo. Shown are an NFW 
halo profile p(r), and a homologous satellite profile a(£) properly shifted to the left and 
upwards. The horizontal dashed arrows refer to deposit based on if>(a) = 1. This would 
steepen the profile, as steep regions of a{£) are deposited at flatter regions of p(r). The 
solid arrows illustrate realistic tidal mass transfer, with a decreasing if) (a). The vertical 
displacements, which grow with r, refer to if) (a) < 1. In this case, the slope at £ is closer or 
even flatter than the slope at r, and more mass tends to be deposited at large r, leading to 
flattening of the halo profile. 

assume that satellite shell £ is being deposited at halo radius r where the slope is a, with 
a given decreasing function if) (a). To visualize this deposit process, we connect by arrows 
two points on the log-log curve a(£) with the corresponding points on the log-log curve p(r), 
illustrating both the (horizontal) distance ratio r/£ and the (vertical) density ratio if). The 
simplified condition if) = 1 (as adopted, e.g., by Syre & White 1998) leads to a correspondence 
via horizontal arrows. In this case, the satellite slope at £ is always steeper than the halo 
slope at r, naturally leading to a steepening of the halo profile. The mean profile evolves 
in this case towards the steepest possible power law, a = 3. On the other hand, when if> 
is properly decreasing with a, the relation between I and r is now illustrated by the solid 
arrows, where the vertical displacement is a growing function of r and a. This is expected 
to weaken the steepening and possibly turn it into flattening of the profile, because the slope 
at £ is now closer to the slope at r, and may be even flatter. Another way to understand 
this flattening effect is by noticing that due to the decrease of if>(a) more satellite mass is 
deposited at outer halo radii. 

We can be a bit more quantitative as follows. With m -C M, we assume that the original 
halo is not affected by tides, and that its mean density profile evolves only by the addition 
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of the satellite mass wherever it is being deposited, namely, 

p fi „ai(r) = p(r) + ^(r)]^, (11) 

where £(r) is defined by mf(r) = m(£). Given that a(r) is minus the logarithmic derivative 
of p(r), it is straightforward to show that the change of a at r due to the merger is 

A , , F(r) dF(r) s a(£) £ 3 

Aar = - J rr^^ 1 , F(r) = -H-—. 12 

Note that the first term in F(r) is l/if>(r), and that F(r) = m(£)/M(r). 

Every power law is a self-similar solution, where Aa(r) = 0. This is because when 
a(r) is constant, so is if>(r), namely a(£)/p(r) is constant, implying that £/r = const., and 
therefore the derivative dF/dr vanishes. However, Aa(r) = does not guarantee a stable 
solution to which the process would converge. 

When a is varying with r, as long as a is small and if) is close to unity, the ratio £/r 
is a decreasing function of r. Thus, if ^[a(r)] =const., the £ 3 /r 3 term determines a positive 
Aa(r), namely continuous steepening towards a — 3. On the other hand, when if) (a) is 
properly decreasing with a, l/ip[a(r)} is an increasing function of r, an increase which can 
balance the decrease of £ 3 /r 3 and produce a fixed point of Act = at some asymptotic slope 
a as . This is a stable solution, where Aa > at a < a as , and Aa < at a > a as . 



4.3. Linear Perturbation Analysis: Summary 

A rigorous perturbation analysis, for a given if) (a) and m/M, is described in the following 
section. We provide here a summary of this analysis and its main result so that the reader 
who is not interested in the mathematical details can skip §5, which is rather technical. 

The linear analysis is based on the following two assumptions, which are argued to be 
of quite general validity: 

• Near the asymptotic solution the slope a(r) is varying slowly, namely, the first and 
higher derivatives of a(r) are much smaller than a(r). This allows us to write a(r) ~ 
cko + eS(r), and examine the behavior of Aa at first order in the small parameter e. 

• The (very small) varying part of a(r) is approximately linear in lnr, namely, the second 
and higher derivatives of a(lnr) are much smaller than da/dlnr. 
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Fig. 6. — The asymptotic slope a as due to a sequence of mergers with a fixed mass ratio m/M. 
The solid curve is the prediction of linear theory, eq. (14), for ACDM (n = —2.33). The 
symbols represent the asymptotic results of the non-linear non-parametric toy simulations 
(§7), demonstrating the success of the linear approximation in the relevant range of m/M. 
The assumed mass-transfer model in both cases is ip{a) = 0.5 /a. The linear results for two 
other power spectra are shown as dashed curves with n marked. The values of a as obtained 
by a distribution of satellite masses [§6, eq. (41)] are marked by the horizontal bars. The 
minimum mass ratio imposed by the dynamical-friction time scale [§6, eq. (39)] is indicated 
by the vertical bar. 

This analysis leads to an explicit expression for Aa at r: 



Here a = a(r), a'(r) = da/dr and ip'{u) = dip/da. Thus, the sign of Aa(r) is the sign of 
the expression in big brackets, and the asymptotic slope a as is the solution of the equation 



The factor in front of the big parentheses in eq. (14) determines the rate of convergence to 
a as ; note that the rate slows down as a approaches a constant and a' — > 0. Recall that ip(ct) 
is given, and that the factors g and A are functions of the mass ratio in the merger and the 
cosmological power spectrum, eq. (10). 




(13) 




(14) 
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The lines in Fig. 6 show the asymptotic slope a as as predicted by the linear theory, 
eq. (14), for a sequence of mergers with the same mass ratio m/M, as a function of m/M. 
The assumed mass-transfer recipe in both cases is if) (a) = 0.5/ a. The middle curve is for 
ACDM, n = —2.33, and the other curves are for n = — 3 and n = —1.67. 

The discussion so far referred to a sequence of mergers with the same mass ratio. In 
§6 we will generalize the analysis to the case of a realistic cosmological distribution of mass 
ratios. 



5. LINEAR ANALYSIS OF THE ASYMPTOTIC PROFILE 

This section is rather technical, and can be skipped by the non-practitioner reader 
without hurting the flow of the discussion. 



5.1. The General Case 

The assumptions and analysis become more transparent if we switch to logarithmic 
variables, 

f = lnr , p = lnp . (15) 

Hereafter we use a tilde to mark logarithmic quantities in general. Then eq. (11) becomes 



Pfinal(^) = p(f) + ln[l + e 



a(i)-p(f)+3(i-f)i 



(16) 



Using the log variables, the local slope is simply (minus) the first derivative of the density 
profile, 

a(r) = — dp(f)/df. (17) 
Then, from eq. (16), the change in the slope in the merger, afi na i(^) — oc{f), is 



Aa(f) 



1 _|_ e p(r)-a(£)+3(r-e) df 
1 



[&(£)- p(r)+3(£-r) 



1 _|_ e p(r)-<r(e)+3(r- 



a(£- A)— -air) -3 — - 1 
dr \ dr 



(18) 



In the second equation we have assumed that the profiles of the halo and satellite are ho- 
mologous, eq. (9). 

In order to make the above expression for Act useful, we need to know the function 
£(f), and especially its derivative d£/dr. This is provided by the mass-transfer prescription, 
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eq. (4), which in the log variables takes the form 

/5(f) - ~a{l) = $(a) , (19) 

with the function ip(a) given, and a = a(f). 

Substituting in eq. (19) the scaling relation for homologous halo and satellite, eq. (9), 
we obtain the equation for mass transfer in the homologous case: 

p(r)-p(£-\) = Q + j>(a). (20) 

Once the scaling factors A and g are given, and the functions ip(a) and a(r) are known, this 
equation should allow us to evaluate the desired 1(f), to be used in eq. (18). 



5.2. First-Order Approximations 

For any pure power law, a(f) = const., we have Act(f) = for all f's. This is because 
then if) is constant, meaning that p—b is constant, so l—f is constant, and then the derivative 
vanishes in eq. (18). However, not all these solutions are stable. In order to find a stable 
solution, we examine Acn(f) for small perturbations of order e about a = cuo, and see how 
the first order in e behaves. 

We first wish to evaluate d£/df at first order in e. Based on the definition of a, the 
homologous mass-transfer condition, eq. (20), can be written as 



I 



ds a(s) = Q + ip(a) . (21) 



Into this equation we substitute the perturbed quantities 

a(f) = a + eS(f) , (22) 

1(f) = £ (r) + ^i(f) , (23) 

where e is small, 5 is of the order of a, and £i is of the order of £. We now collect terms of 
the same order. The zero's order term yields 

[e (r)-\-r\a = Q + $(a), (24) 

h{f) = [Q + i>(a )]/a +X + f . (25) 

The first order term gives 

t0 \s5(s)+£ 1 (f)a = ^^5(f), (26) 
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ii(f) = — 

a 



5(f) 



e -x 



dsS(s) 



(27) 



where ip'(a) = dip /da. Using these results we can now calculate d£(r)/dr to first order: 



dr 



e 

OiQ 



ipl m+S[f) . s{ i a -\) 



(28) 



where 5'(r) = d<5/dlnr. 

When we plug this result into eq. (18), we obtain after some algebra 

Aa(f) = 



1 _|_ e p(f)-a(e )+3(f-e ) 
\ip{a ) a 



(29) 



x 



+ 0{e 2 ) . 



_ip(a ) 

We can further simplify the above expression by substituting £q with its explicit value taken 
from eq. (25). Using eq. (20) and the definition of A, we obtain 

e 



Aa(f) 



1 + \- 3 g- 3 / a o^(a ) 1 ' 3/a ° 



x 



i/j(a ) 



8'(f) + 5(f) -5(f+[g + ^(a )]/a ) 



(30) 
+ 0(t 2 ) . 



Since the above expression is first-order in e, we may safely replace ao by a(f) because 
the difference between the two is also first order, eq. (22). Additionally, eq. (22) implies the 
identities 

e( J'(f) = a'(f) , (31) 
e5(f) - eS (r + [g + ijj(a )]/a S j = a(f) - a (f + [g + ^(a )]/aio) , (32) 

which allow us to rewrite eq. (30) in terms of a(f) and its first derivative a' if) [with no other 
explicit functions of f such as 5(f) or 5'(f)]: 

a'(f) 



Aa(f) 



1 + \- 3 g- 3 l a i)(a) l - 3 l a 



(33) 



ip(a) a I tp(a) 



a [f +[g + ip(a)]/aj — a 
a'(f) 



+ 0(e 2 ) . 



Here, and in what follows, a stands for a(r). This notation hides the explicit dependence of 
Aa(f) on e, but for the calculation to be self consistent e has to be on the order of 



a 



e ~ 



(f+[g + i>(a)]/o^ 



a 



a 



(34) 
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We thus require the right-hand side of eq. (34) to be much smaller than unity. 

In principle, the sign of Aa should allow an analysis of the evolution of a into a fixed 
point. However, the third term in eq. (33) is still problematic; it involves the function cc(-) 
evaluated recursively at a point which depends on a itself. In order to make this term useful, 
one must make further assumptions regarding the functional form of «(•). 



5.3. The Linear Regime 

The most natural assumption for the functional form of a(f) is that it is linear in f. 
As seen in Fig. 3 of DDH, this is clearly a good approximation for the profile used in our 
merger N-body simulation, and it is also valid for the family of profiles described in eq. (1). 
In other words, we assume that the second and higher derivatives of a(r) are much smaller 
than the first derivative — at least in the interval [f, — A]. Under this assumption we may 
approximate the above problematic term by 



a + [g + ip(a)]/a^ - < \ ( r ; 



1 



a'(f) a 
which brings eq. (33) to finally become 

a'(f) 



-[Q + ip(a)\, (35) 



Aa(r) 



1 + \- 3 Q-y a ?p(a) 1 - 3 / a 

(a(a- 3)4£r + 3[£ + ^(a)] 1 . (36) 



In this approximation, the sign of Aa(r) is the sign of the expression in big brackets, 
and the asymptotic fixed point slope a as is obtained by solving an explicit equation of a: 

a{a - 3)^1 + 3[g + i>{a)] =0. (37) 

The factor in front of the big brackets of eq. (36) affects the rate at which ct(f) converges to 
the asymptotic slope. This concludes the derivation of eq. (14) and eq. (14) introduced in 
§4.3 for the case of a given mass ratio m/M. 



6. A COSMOLOGICAL DISTRIBUTION OF SATELLITES 



The analysis so far referred to mergers with a fixed mass ratio m. A real sequence 
of mergers consists of a cosmological mix of satellite masses. The probability for a mass 
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ratio m in a merger can be estimated using the Extended Press Schechter scheme (Lacey & 
Cole 1993, EPS), which has been tested and calibrated using N-body simulations (Lacey & 
Cole 1994; Somerville et al. 1999, Fig. 2). We have performed Monte Carlo realizations of 
merger trees based on EPS (using a scheme developed by Yuval Birnboim, to be published 
elsewhere, based on the scheme of Somerville & Kolatt 1999), and found that the probability 
distribution of m is quite insensitive to the actual power spectrum of initial fluctuations, and 
it has no obvious time dependence. The fraction of number of mergers with mass ratio m < 1 
in the interval (m, m + dm) can be approximated by 

p(m) dm oc m " 3 / 2 dm. (38) 

The relevant distribution of mass ratio for the buildup of the inner profile should also 
take into account the time it takes each satellite's orbit to decay by dynamical friction from 
the halo virial radius R to the central region. The duration of the decay process can be 
estimated for circular orbits in an isothermal halo (as in BT, eq. 7-26). We obtain 

t 1.17 R 0.18 

mm(l/ra) V circ m\a{l/m) 

where t{ is the Hubble time when the merger starts with the satellite at the halo virial radius, 
and where an Einstein-deSitter cosmology is assumed at that time. On one hand eq. (39) is an 
underestimate because we have ignored the weakening effect of mass loss on the dynamical 
friction. If we adopt the mass-loss recipe m(r) oc M(r) (which can be deduced from the 
resonance condition for isothermal halo and satellite) we obtain that tf r - lc becomes about 
twice as long as when mass loss is ignored. On the other hand eq. (39) is an overestimate in 
the case of eccentric orbits, where dynamical friction is strong already at early stages of the 
merger. In our merger simulations we find that a satellite on a typical orbit of eccentricity 
1:6 decays to the center on a time scale shorter by a factor of ~ 3 compared to the decay 
time of a similar satellite on a circular orbit. Thus, the two effects roughly balance each 
other, and we can adopt eq. (39) as a crude approximation. Note that for m ^ 1 the estimate 
in eq. (39) breaks down (e.g., the log term drives tf ric to grow with m for m > e _1 ). Since 
we wish to push the crude analysis to include such mergers despite the breakdown of our 
approximations there, we simply keep tf r i c constant in this range. 

The satellite is irrelevant for the cusp at to if ^fric(^) > to — ti, which therefore defines a 
minimum relative satellite mass m min . For larger m we crudely multiply the mass function 
of eq. (38) by the approximate correction factor /f r i c (m) = [1 — t {ric (m) / (t — t{)]. For an 
order-of-magnitude estimate, we assume that a typical merger occurs at U — to/2 (e.g., for 
an isothermal halo this is the typical time of collapse of the inner half mass of the halo). 
This implies m min ~ 0.067 and f ilic ~ (1 — 0.18/[raln(l/m)]) for m min < m < e~ l . A lower 
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Fig. 7. — The asymptotic slope due to a sequence of mergers with a cosmological distribution 
of mass ratios, as obtained in the linear approximation by solving eq. (41), shown as a 
function of the fluctuation power index n. 

value for if ric , or an earlier ti/t , would correspond to lowering the factor 0.18 in f {ric and 
therefore to a smaller m m i n . Since m m ; n ~ 0.067 is already much smaller than unity, the 
effect of the actual value of m m j n on our final result is weak. 

Given the corrected probability distribution /f ric (m)p(m), we compute the average over 
m of Aa(r, m) as given in eq. (14). We obtain the asymptotic slope a as for the distribution 
of satellites by requiring that the average increment vanishes: 



Note that we have expressed g and A in terms of m and v, using the scaling eq. (10). The 
factor a'{f)/a 2 has been taken outside the integral and dropped from the equation because 
it is independent of the stochastic variable m of the current merger (though it is affected by 
the mass ratios in previous mergers throughout the halo history). 

Fig. 7 shows the value of the asymptotic profile obtained by solving eq. (41) for the 
above distribution of mass ratios and with ip(a) = 0.5/ a as a function of the power index n 
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of the fluctuation power spectrum [v = (3 + n)/2]. Two values are also marked as horizontal 
bars in Fig. 6. For ACDM, with n ~ —2.33, we obtain a as = 1.29. As seen in Fig. 6, 
the value of « as for a distribution of mass ratios is similar to the values obtained from 
eq. (14) for to ~ 0.32 [which turns out to be close the average (to) ~ 0.33 of the distribution 
/fric(^)p( r ^)] ) indicating that a typical merger in terms of its effect on the cusp profile is 
with to ~ 1/3. This is true for any n. 

The reason for the robustness of the dominance by to ~ (to) ~ 1/3 mergers and the 
insensitivity to the actual value chosen for m m i n is the relative flatness as a function of to of 
the term multiplying the expression in big parenthesis in the integrand, eq. (41). The factor 
/f r i c (m) is roughly constant in most of the range, while p(m) oc m~ 3 / 2 . The term expressed 
as a big fraction is roughly proportional to m 1 ^ 1-3 /") which is about oc m 5 to to in the 
relevant range of v and a. Together they introduce a moderately decreasing dependence on 
to, which does not become too large even at small to, as long as m min is not significantly 
smaller than 0.01, say. 

In Fig. 7 we see a gradual variation of a as with n: for the extreme case of n = — 3 we 
get a minimum slope of a as = 0.98 and for n = —1.67 we obtain a as = 1.68. This range is 
comparable to the variation seen in simulations of a single cosmology, e.g., ACDM. We shall 
see below (§7) that the measured cusp slopes are actually expected to be somewhat less n- 
sensitive than a as . A systematic trend of such magnitude is hard to measure in cosmo logical 
simulations of different power spectra. There are hints for a similar trend already in the 
simulations of Cole & Lacey (1996, Fig. 9), and SCO report from their simulations an 
average of a m = 1.3 ±0.07 for n = —2, then a- m = 1.6 ±0.09 for n — — 1 and a in = 1.8 ±0.09 
for n — 0. 

As a side complementary argument, we use a very simplistic toy model to derive a crude 
upper limit for a m . Following SCO, we consider the accumulation of undigested satellites in 
the halo center. Based on the assertion that the satellites' inner densities are significantly 
higher than that of the host halo, one assumes in this model that the satellites accumulating 
in the halo center keep their original profiles unchanged by tidal effects. The mean halo 
profile that is built by a cosmological distribution of satellites is then 

p(r) = / dmf ilic (m)p(m)a(m;r) (42) 

where <t(to; r) is the profile a(£) at r = £ of a satellite of mass ratio to, and where the to 
dependence is given by the cosmological scaling for the given power spectrum. This integral 
is easy to evaluate if we assume that /f r i c (m) = const., adopt p(m) oc m~ 3 / 2 as in eq. (38) 
based on the realizations of the EPS model, and crudely take a(m; r) to be a step function 
of height <7 S out to a radius £ s , both scaled with to as in eq. (10). At a given r, there is 
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a in =£* in +Aa(r=0.01) 

f=0.b/a 

m/M=l/3 




20 40 60 80 100 

merger 

Fig. 8. — Evolution of the cusp slope in a sequence of mergers with homologous satellites of 
m/M = 0.33. In this version of the toy simulation the functional form of eq. (1) is enforced 
with a out = 3, starting with either a in = 1 or 2. The change in a in is assumed to be given by 
A«(r) at r = 0.01r s . The cusp slope converges slowly to an asymptotic value, as predicted 
by the linear analysis. 

a contribution to p(r) from all satellites of £ s > r, which translates by eq. (10) to a lower 
bound for the integration at a minimum mass m = r 3 ^ 1+u \ The result, with v — (3 + n)/2, 
is 

a- m = 3(4 + n)/(5 + n). (43) 

SCO argued for a flatter slope, a- in = 3(3+n)/(5+n), which could have resulted from eq. (42) 
if one had assumed fhi C {m) p(m) oc m~ x . Based on the above discussion of dynamical friction 
we expect /f r i C (m) to be a moderately increasing function of m for m > m min . This would 
flatten the slope predicted in eq. (43) and thus make it an upper bound for the actual cusp 
slope. An effective behavior of /fric(^) oc m 1//3 , with p(m) oc m _3//2 , would roughly recover 
the n dependence found in the simulations of SCO. This simplified toy model ignores the 
additional effects of tidal mass transfer [which we model in the current paper via ip{ a )} an d 
the associated effects of tidal distortion of the satellites' inner regions, which may go either 
way. Nevertheless, we notice that eq. (43) indeed overestimates our model prediction for a as 
for any n, and thus serves as an upper bound to a- in . 
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7. TOY SIMULATIONS 

In order to test the predictions of the linear analysis of §5 for the asymptotic profile, 
and to study the evolution towards this asymptotic profile, we perform toy simulations of 
the profile buildup by a sequence of mergers, where we implement straightforwardly the 
mass-transfer recipe of eq. (4) with vp(a) = 0.5/ a in the range a > l. 3 We start with a 
halo of a given mean density profile p{r) and the corresponding monotonically increasing 
a(r). We consider a merging satellite of mass ratio m/M with a homologous profile properly 
scaled using eq. (10). To materialize the merger we solve eq. (4) numerically for £(r), with 
the possible modification for deposit without stripping at a < 1. We then assume that the 
halo profile changes only due to the added stripped satellite material, eq. (11), and compute 
the exact change in slope Aa(r) at any r using eq. (12). The following merger is performed 
with a satellite that is a scaled version of the new halo profile, and so on. 

As a first simple test, we enforce at all times a mean density profile obeying the functional 
form of eq. (1), with r s and a ou t fixed at 1 and 3 respectively. The inner slope a in starts at an 
arbitrary value in the range 1 < a m < 2 and is allowed to change as a result of the merger. 
We assume here that a m changes by Aa(r) as evaluated at some small radius r. Fig. 8 shows 
the resultant evolution of a m for m/M = 0.33 and r = 0.01r s , when starting alternatively 
from either a in = 1 or a m = 2, and for three different power spectra of fluctuations. We see 
that, indeed, there is a slow convergence to an asymptotic profile. The obtained asymptotic 
slopes are close to the predictions of the linear analysis, Fig. 6. The slight deviations of less 
than 10% are mostly due to the fact that we compute Aa at a finite small r rather than at 
r -> 0. 

We next perform a non-parametric toy simulation where we do not enforce a specific 
functional form and do not constrain the outer slope to remain steep. We only make sure 
that a(r) is monotonic and relatively smooth, as required for the implementation of the 
deposit scheme. At each step, the mean density profile of the halo, p(r), is stored in an array 
of shells spaced logarithmically (1000 shells, Alnr = 0.05), and spline interpolated into any 
desired r. The slope a(r) is computed via 4'th order spline interpolation. To make sure 
that a(r) remains monotonic and smooth, a(r) is first smoothed with a Gaussian of width 
A lnr = 0.5, and then interpolated inside intervals of A lnr = 2 (with an overlap of ±1 with 
the neighboring intervals) using locally a functional form for a(r) corresponding to eq. (1) 
with all the parameters free. This is a general smoothing procedure that is not sensitive to 



3 Certain modifications to the recipe may be needed near and below a ~ 1, where our adopted ip{ a ) i s 
not well established, and where deposit may occur without stripping. For example, if eq. (4) yields a deposit 
radius r < £, one may wish to replace it by r — £, assuming that the surviving satellite supports itself in a 
finite configuration while its center sinks to the halo center. 
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Fig. 9. — Non-parametric and unconstrained toy simulation evolution of slope profile a(r) 
due to a sequence of mergers (n = 1,400) with mass ratio m/M = 0.3. The initial profile 
is eq. (1) with a- m = 0.5 and a ut = 3. A power-law region develops below the radius where 
Aa = 0, with a slope that grows slowly from near unity to the asymptotic value. For a long 
time, the profile maintains a shape of an inner cusp and a steeper outer halo, as seen in 
cosmological simulations. 



the use of this specific functional form in the local fit. Our initial profile is again given by 
eq. (1), with r s put at the center of the sampled lnr range. 

The ultimate test of the linear approximation is shown in Fig. 6, where the results 
of the non-parametric toy simulations are presented by the symbols on top of the curves 
describing the linear predictions. We see that the actual sequence of mergers with a fixed 
mass ratio m/M converges to an asymptotic slope very similar to that predicted by the 
linear approximation, for all values of m/M in the range of interest. Significant deviations 
are obtained only for very small m/M values, below any realistic value of m min as implied 
by dynamical friction time-scales. 

Fig. 9 shows how a(r) actually evolves in the non-parametric toy simulations where 
no constraints are applied on the outer slope. When a in is below unity there is a rapid 
evolution towards a ~ 1 in the range r < r s . This is partly due to the factor 1/a 2 driving 
Aa in eq. (12) and partly due to the processes explained in DDH as reflected in the shape 
adopted for ip(a) at a < 1. A power-law region with a ~ 1 develops below the radius where 
Aa = 0, and it quickly extends to smaller radii. The value of a in this inner power-law 
region grows slowly from near unity towards the asymptotic value a as . In parallel, since 
no constraints are applied on the outer slope, it can decrease slowly from a out = 3 towards 
a as . The shape of the profile during the long interim period resembles a cusp-halo structure 
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qualitatively similar to the NFW-like profiles of halos seen in cosmological simulations. It 
takes 10 to 20 mergers (of m/M = 0.33), during which the halo mass grows by about two 
orders of magnitude, for the shape of the profile to strongly deviate from this characteristic 
shape. In reality, we may expect the outer profile to keep a steep slope in response to other 
physical processes not modeled here (see §1), and therefore the asymptotic profile to resemble 
a shape qualitatively similar to eq. (1) with a m = a as . 

8. DISCUSSION AND CONCLUSIONS 

We have addressed the origin of a cusp in dark-matter halos within the hierarchical 
clustering scenario, based on tidal effects during the halo buildup by merging of relatively 
compact satellite halos. In Dekel, Devor & Hetzroni (2003) we highlight the rapid steepening 
of a flatter core of slope a < 1 into a cusp of a > 1 as a result of vanishing tidal mass transfer 
in such a core. Merger N-body simulations indicate here that the cusp is quite stable under 
a sequence of such mergers, we have developed a simple prescription for tidal mass transfer 
in mergers, argued using N-body simulations that it provides a useful approximation, and 
showed that it naturally leads to a stable cusp with an asymptotic profile of slope a as slightly 
larger than unity. Our toy model is one way to explain the cusps seen in cosmological N-body 
simulations of the hierarchical clustering scenario, and it provides a tool for addressing other 
processes which may explain the observed flat core in some galaxies. 

The mass-transfer recipe relates each radius I in the initial satellite with a deposit radius 
r in the halo via the ratio of mean densities p(r)/a(£) = ip[a(r)]. We argued based on quali- 
tative considerations and found in merger N-body simulations that ip(a) has a characteristic 
shape which is not too sensitive to the details of the merger orbit; ip(a) is a decreasing 
function of a and it obtains values significantly smaller than unity. A crude fit seems to be 
provided by ip(a) = 0.5/a in the range 1 < a < 3. This makes the mass transfer more effi- 
cient than what might be inferred from the naive resonance condition, ip = 1, and gradually 
more efficient where the halo profile is steeper. Beyond our immediate purpose of studying 
the formation of a cusp, this mass-transfer recipe provides a useful tool for general studies 
using semi-analytic models of galaxy formation. 

Admittedly, we have only tested in this paper the robustness of ip(a) to the merger orbit 
type. The choice of typical properties for the simulated halo and satellite, which relate to 
each other in general agreement with the scaling relation expected in the ACDM cosmology, 
makes the obtained ip(a) a sensible crude guess. Nevertheless, the robustness of ip(a) to 
variations about the simulated properties of halo and satellite must be tested using a more 
comprehensive set of simulations. We have adopted ip(a) = 0.5/a as an illustrative example 
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for obtaining actual numerical values for a as , but the convergence analysis, which is the heart 
of this paper, is much more general. 

After understanding why a decreasing if) (a) should lead to a stable asymptotic profile, 
we have performed a linear perturbation analysis of this process, and confirmed its validity 
by toy simulations of tandem mergers based on the adopted mass-transfer recipe. The 
asymptotic slope a as is typically somewhat larger than unity. The system develops a cusp 
of a ~ 1 which slowly grows towards a as while the outer halo maintains a steeper profile 
for a long time. For a typical cosmological distribution of merging satellites in the ACDM 
cosmology (n = —2.33) the asymptotic slope is a as ~ 1.3. This is close to the typical cusp 
slope found in cosmological N-body simulations (e.g. Power et al. 2002 and other references 
in §1). The typical mergers driving this process are with a relatively high mass ratio of 
m/M ~ 1/3. This means that one should not take our quantitative results too literally 
because our toy model assumes m <C M. 

The asymptotic profile obtained in our model has a relatively weak dependence on the 
power spectrum of fluctuations. In particular, the effect of tidal compression (§2) and the 
1/a 2 factor in Aa [eq. (12)] drive the profile to a ^ 1 independently of the fluctuation power 
spectrum. Also, the cosmological dependence of the ultimate asymptotic slope is weakened 
by the fact that the profile is determined by mergers of relatively large mass ratio. For values 
of power-spectrum index in the range —3 < n < — 1 the asymptotic slope is predicted to 
vary gradually in the range 1 < a as < 2. This can be compared to the only slightly weaker 
trend indicated in the cosmological N-body simulations of SCO. 

Syer & White (1998) also addressed the profile resulting from mergers and obtained 
results somewhat different from ours. They adopted the condition of tidal stripping at 
resonance, if) — 1, ignoring the a dependence of if). When we substitute the recipe if) — 1 
in eq. (12), we find that Aa is always positive, for any merger and at any r. Indeed, when 
trying to repeat the Syer & White toy simulations with higher resolution and following more 
mergers we actually find that the profile does not really converge to a stable cusp but rather 
continue to steepen slowly towards a- in = 3. Only when using the revised mass-transfer 
prescription where if) is decreasing with a do we obtain convergence to a flatter asymptotic 
profile. 

Syer & White (1998), Nusser & Sheth (1999) and Subramanian, Cen & Ostriker (2000) 
describe different toy models which predict a m = (9+3n) / (5+n), 4 recovering the hierarchical- 
clustering toy model by Peebles (1980, eq. 26.8). While this scaling relation may be in 
reasonable agreement with the simulations (and with our model predictions) for n > — 1, 



4 Nusser & Sheth also predict an upper bound of a; n = (9 + 3n)/(4 + n). 
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it admits values below a m = 1 for n < — 2, in conflict with the findings of the simulations. 
Our revision of the toy model of SCO yields a in < 3(4 + n)/(5 + n), to be flattened further 
by the m dependence of dynamical friction. The basic model analyzed in the current paper 
does not admit cores significantly flatter than a m ~ 1 (that is by non-dissipative processes 
only) and is thus in better agreement with the simulation results. 

The gravitational processes leading to a cusp are modeled in our analysis as tidal effects 
during the buildup of the halo by a sequence of mergers. This picture is likely to be valid 
in the CDM hierarchical clustering scenario, where numerous sub-galactic halos exist on all 
scales and are continuously merging (e.g., Klypin et al. 1999b; Moore et al. 1999a; Springel 
et al. 2001). This merger picture is confirmed by a careful inspection of high resolution CDM 
simulations, where a special effort is made to identify the merging clumps which otherwise 
could have been easily missed (e.g., Wechsler, Dekel et al. , in preparation). It is therefore 
clear that cores, independently of how they form, cannot survive in a pure CDM scenario 
without some modification — they efficiently turn into cusps as a result of the mergers with 
the relatively compact building blocks. 

Nevertheless, a cusp, though somewhat flatter, is reported to be seen also in simulations 
where the initial fluctuations had less power on small scales, thus suppressing the number 
of sub-galactic satellites and the associated merger rate (Moore et al. 1999b; Avila-Reese et 
al. 2001; Bullock, Kravtsov & Colin 2002). Further indications for the generality of cusp 
formation comes from simulations by Huss, Jain & Steinmetz (1999), who find that cusps 
also form as a result of collapse from spherical initial density perturbations with suppressed 
random velocity perturbations, as well as from the simulations of Alvarez, Shapiro & Martel 
(2002), who find that cusps arise from the gravitational instability and fragmentation of 
cosmological pancakes. One way to explain this is by noticing that the asymptotic cusp 
formed in the CDM scenario is driven by mergers with relatively massive satellites (of typical 
mass ratio 1:3, see Dekel et al. 2002), and realizing that such mergers do happen even when 
small-scale power is suppressed and when pancakes fragment. It is therefore possible that 
the cusp is actually driven by mergers to a certain extent also in some of these cases. Our 
toy model based on mergers would not be directly applicable in cases where a cusp arises 
while the halo forms by smooth accretion. The gravitational processes involved in such cases 
somehow mimic a behavior similar to the merger case. We note, for example, that the tidal 
compression discussed in DDH is expected to amplify density perturbations and possibly 
make them behave in certain ways like merging clumps. 

However, one should accept in general that the halo buildup is a complex gravitational 
process, whose different aspects can probably be modeled in more than one way, e.g., as a 
violent relaxation process driven by fluctuations in comparison with a sequence of mergers 
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and substructure accretion. 5 There is also the intriguing possibility that general statistical 
considerations in phase space may provide some clue for the origin of a universal halo profile 
(e.g., Taylor & Navarro 2001). Unfortunately, we currently know of no viable alternative to 
the tidal effects in mergers as a simple model for the origin of a ^ 1 cusps in the cosmological 
simulations. The idealized merger picture provides one possible toy model, approximately 
valid at least in hierarchical clustering scenarios, within which we understand the origin of 
the cusp in simple terms. However, the complexity of the problem suggests that this is not 
the final word on the issue. 

Our result implies a necessary condition for the survival of cores in halos independent of 
their origin — that satellites should be prevented from adding mass to the halo cores. This 
could be avoided in a CDM scenario if feedback processes manage to puff up the satellites 
and make them disrupt before they merge with the halo cores (as proposed in DDH). We 
have not explicitly addressed in this paper the sufficient conditions for the formation of 
cores, but one can imagine that a sequence of mergers with low-density satellites, where 
the mass is predominantly deposited outside the inner halo region, would indeed flatten the 
inner profiles (work in progress). Other processes may also contribute to the development of 
halo cores. The proposed scenarios include, for example, the disruption of cusps in merging 
satellites by massive black holes (Merritt & Cruz 2001), the heating by gas clouds spiraling 
in due to dynamical friction (El-Zant, Shlosman & Hoffman 2002), the angular-momentum 
transfer from a big temporary rotating bar (Weinberg & Katz 2002), and the delicate resonant 
reaction of halo-core orbits to the tidal perturbation by the satellite, which could be a strong 
effect if the dark-matter distribution is much smoother than the current state-of-the-art N- 
body simulations (Katz & Weinberg 2002). We stress again that, no matter what the origin 
of the core might be, our analysis implies that such cores could survive only if they are not 
perturbed by significant mass transfer from merging compact satellites, which implies that 
small CDM halos must be puffed up before they merge into bigger halos. We also stress that 
none of these scenarios seem to be capable of explaining the formation of cores in halos on 
the scales of clusters of galaxies within the CDM scenario. 

Supernova feedback effects are probably not strong enough for turning cusps into cores 
in halos with rotation velocities higher than ~ 100 km s" 1 (Geyer & Burkert 2001; Gnedin & 
Zhao 2002). However, they are possibly sufficient for the necessary indirect puffmg-up of the 
merging satellites. Gnedin & Zhao (2002) estimate that direct feedback effects may reduce 
the central satellite densities by a factor of 2 to 6. Based on our simulations with puffed-up 



5 We encounter an analogous duality, for example, when one manages to explain the origin of galactic spin 
alternatively via tidal torque theory applied to shells and by summing up the orbital angular momenta in a 
cosmological sequence of mergers (see Mailer, Dekel & Somervillc 2002). 
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satellites, this may be enough by itself to avoid the steepening from a core to a cusp in 
the framework of the merging scenario. However, if the (so far inconclusive) observational 
clues for cores in clusters of galaxies are confirmed, supernova feedback is unlikely to provide 
a viable explanation for their origin. In this case one may search for higher efficiency in 
supernova feedback either due to microscopic effects such as porosity in a muliphase ISM 
or due to hypernova from very massive stars (Silk 2002). Alternatively, one may appeal to 
stronger feedback mechanisms, perhaps associated with radio jets from AGNs. This process 
may be indicated by an observed correlation between AGN activity and bright galaxies 
in SDSS (Kauffmann et al. 2003, in preparation), and it may be needed independently 
in order to explain the missing baryons in big galaxies (Klypin, Zhao & Somerville 2002) 
and in clusters. Otherwise, such large cores may present a real challenge for the standard 
hierarchical clustering scenario 

The cusp/core problem is only one of the difficulties facing galaxy formation theory 
within the CDM cosmology. It turns out that other main problems can also be modeled 
by tidal effects in mergers, and may also be resolved by the inevitable feedback processes. 
For example, Mailer & Dekel (2002) addressed the angular-momentum problem, where sim- 
ulations including gas produce disks smaller than the galactic disks observed (Navarro & 
Steinmetz 2000 and references therein; Governato et al. 2002), and with a different internal 
distribution of angular momentum (Bullock et al. 2001b; van den Bosch, Burkert & Swaters 
2001). A toy model has been constructed for the angular-momentum buildup by mergers 
based on tidal stripping and dynamical friction, which helps us understand the origin of 
the spin problem as a result of over-cooling in satellites. A simple model of feedback has 
then been incorporated, motivated by Dekel & Silk (1986). This model can remedy the 
discrepancies, and in particular explain simultaneously the low baryon fraction and angular- 
momentum profiles in dwarf disk galaxies. 

Various feedback effects may also provide the cure to the missing dwarf problem, where 
the predicted large number of dwarf halos in CDM can possibly match the observed number 
of dwarf galaxies only if the mass-to-light ratio in these objects is very high (Klypin et 
al. 1999b; Moore et al. 1999a; Springel et al. 2001; Kochanek 2001). Bullock, Kravtsov 
& Weinberg (2000), Somerville (2002) and Tully et al. (2002) appeal to radiative feedback 
effects which prevent the formation of small dwarfs after cosmological reionization at z ~ 7, 
Scannapieco, Ferrara & Broadhurst (2000) and Scannapieco & Broadhurst (2001) address the 
destructive effect of outflows from one galaxy on neighboring protogalaxies via ram pressure, 
and Dekel & Woo (2002) study the role of supernova feedback in determining the relevant 
global properties of dwarfs and larger low-surface-brightness galaxies. We note that while 
the requirements from feedback in explaining the dwarf-galaxy properties and the angular- 
momentum problem are not too demanding, the solution to the core problem requires that 
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the dark-matter distribution be affected by feedback, which is a non-trivial requirement. 

Nevertheless, the successes of such toy models in matching several independent observa- 
tions indicate that they indeed capture the relevant basic elements of the complex processes 
involved, and in particular that feedback effects may indeed provide the cure to some or all 
the main problems of galaxy formation theory within the ACDM cosmology that does so 
well on larger scales. The alternative solution involving Warm Dark Matter (e.g., Hogan 
& Dalcanton 2000; Avila-Reese et al. 2001; Bode, Ostriker & Turok 2001) seems to still 
suffer from the cusp problem, it may still fail to reproduce the angular-momentum profile 
in galaxies (Bullock, Kravtsov & Colin 2002), and it may be an overkill where the forma- 
tion of dwarf galaxies is totally suppressed once the inevitable feedback effects are included 
(Bullock 2001). The speculative alternatives involving self-interacting dark matter are even 
more problematic (Spergel & Steinhardt 2000; Dave et al. 2001; Hennawi & Ostriker 2002). 
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